Inferring the time-dependent complex Ginzburg-Landau equation from modulus data 



Rotha P. YuQ David M. PaganinQ and Michael J. MorgarH 
School of Physics and Materials Engineering, Monash University, Victoria 3800, Australia 

(Dated: February 6, 2008) 

We present a formalism for inferring the equation of evolution of a complex wave field that 
is known to obey an otherwise unspecified (2+l)-dimensional time-dependent complex Ginzburg- 
Landau equation, given field moduli over three closely-spaced planes. The phase of the complex 
wave field is retrieved via a non-interferometric method, and all terms in the equation of evolution 
are determined using only the magnitude of the complex wave field. The formalism is tested using 
simulated data for a generalized nonlinear system with a single-component complex wave field. The 
method can be generalized to multi-component complex fields. 
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I. INTRODUCTION 

In quantum mechanics only the probability density of 
the system can be directly measured. Measurement of 
the probability density alone provides restricted infor- 
mation about the physical phenomenon under investiga- 
tions. To make sense of the data collected by a measur- 
ing instrument, and hence gain insight into the dynam- 
ical behaviour of the system, it may be postulated that 
the system must obey an, otherwise unknown, physical 
law. These laws of physics are commonly written in the 
form of partial differential equations. To infer the pre- 
cise form of the partial differential equation involves using 
the data to impose constraints on the equation leading to 
the identification of the equation of evolution This 
approach has been limited to real scalar fields or known 
complex scalar fields (see e.g., 0,0>@)- In such cases the 
equation of evolution of a complex scalar field can only 
be identified if both the probability density (intensity) 
and phase of the complex wave field are known a priori. 
However, in general, the phase of the complex wave field 
is not directly measurable, as only the probability den- 
sity of the system can be measured. In such a case, we 
ask a fundamental question: can we infer the equations 
of evolution of a complex system given only the modulus 
data of its complex wave field? 

There has been much recent work directed towards 
non-interferometric retrieval of the phase of a complex 
wave field from modulus data 0, 0, El • However, such 
an approach assumes a known equation of evolution for 
the wave field. Here, we infer the evolution equation by 
simultaneously retrieving the phase of the wave field and 
all terms in a generic evolution equation of the wave field. 
In a previous paper [T^j, we demonstrated the feasibil- 
ity of determining the evolution equations of a complex 
field, given modulus data alone. The work of that paper 
comprised a means for "measuring" a dissipative (2+1)- 



dimensional nonlinear Schrodinger equation. However, 
the approach was limited to a smaller class of systems 
than that treated here. 

One of the most studied equations in mathematical 
physics is the time -dep endent complex Ginzburg-Landau 
equation (see e.g., pd|). In this paper we provide a means 
for "measuring" or inferring the time-dependent complex 
Ginzburg-Landau equation, from modulus information 
(or the probability density distribution of the system) . 

The organization of the paper is as follows. In Sec. ITT1 
we introduce the time-dependent complex Ginzburg- 
Landau equation. Section IIIII discusses wavefunction 
phase retrieval in the presence of diffusion. Section Hvl 
discusses the sensitivity of the retrieved phase on the 
error in the diffusion parameter. In Sec. we ap- 
ply our formalism for extracting all terms in the time- 
dependent complex Ginzburg-Landau equation from sim- 
ulated modulus data. Section IVII generalizes our discus- 
sion to multi-component complex fields, and in Sec. IVIII 
we conclude with a discussion of some implications of this 
formalism and suggest future directions. 



II. THE TIME-DEPENDENT COMPLEX 
GINZBURG-LANDAU EQUATION 

The time-dependent complex Ginzburg-Landau (TD- 
CGL) equation has the form 



ia£ + (1 - w?)Vi + /(I) + ig{I) 
oz 



* = 0, (1) 
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where VP = ^(x,y,z) is the complex wavefunction, I = 
I(x, y, z) — \$!{x, y, z)\ 2 is the probability density, or the 
intensity, a and rj are real parameters, /(/) and g(I) 
are generalized nonlinear real- valued functions of /, z is 
the evolution parameter and Vj_ = (d/dx,d/dy) is the 
gradient operator in the x-y plane. Special cases of the 
TDCGL equation describe a wide variety of both clas- 
sical and quantum systems, such as monoenergetic elec- 
tron beams [12j , beamlike monochromatic scalar electro- 
magnetic waves ^3 ' intense scalar electromagnetic fields 
in nonlinear media |14|. Bose- Einstein condensates [l5| . 
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uncharged superfluids, and vortices and strings in field 
theory |3. 

When we refer to "measuring" or inferring the equation 
of evolution of the TDCGL equation (|TJ, we shall mean 
solving for all the parameters (a and rf), functions (/(/) 
and g(I)) and the phase arg('i') of the complex scalar field 

from given modulus data information /. In measuring 
the TDCGL equation (JTJ , we limit ourselves to simulated 
modulus data. Consequently our work proceeds in two 
stages. In the first stage we focus on generating simu- 
lated data from the known TDCGL equation (i.e., 
forward evolution of Eq. with a, 77, /(J), g(I) and 
arg('J') specified). The second stage solves for the TD- 
CGL equation based on the assumption that only the 
modulus information I is known. 

In the forward evolution, we solve Eq. Q using 
a fourth-order Runge-Kutta differentiation scheme (see 
e -g- ; [9l fToL IT^|L The simulations were performed with 
a z-stcp Az = 10~ 7 running for 300 iterations. The 
grid size is 1025 x 1025 within the simulation domain 
[0,1] x [0,1], corresponding to the spatial step 1/1024. 
For concreteness, our simulations use a — 1,?/ = 2, 
/(/) = 100sin(7r/) with the power dissipation function 
specified by g{I) = 3J 2 ; however, other parameters and 
functions are of course possible. 

The initial condition is based on a Gaussian intensity 
profile with sinusoidal modulation 



A 



1 + Se »( »'°) cos (2nn(x - x )) 
l + ^e - ^ - ^ 11 ) sin (2irn(y — y )) 
e 2\ w ) 



where A is a constant denoting the peak of the Gaussian 
profile of width W located at r = \J x 2 + y 2 = tq. We 
choose A = 10, W = 8, S = 0.01, r = 0.5 and n = 
20. The initial phase of the system is specified by the 
Gaussian profile 



<f>(x,y) = A$e 



2 V w ) 



(3) 



where ^4$ is a constant, which can be used to vary the 
average phase gradient of the system. We performed sim- 
ulations with A$ in the range [0,2]. It turns out that 
within this range, the average phase gradient (|V$|) is 
also approximately in the range [0, 2]. 

As \& evolves we output the intensity profile to a file at 
every 100 iterations. The profile at the 100th iteration is 
denoted by zq; the 200th iteration is denoted by Z2 and 
the 300th iteration is denoted by Z4 (see Fig. QJ. The 
phase profile of \& is also written to a file at every 100 
iterations for later comparison with the retrieved phases. 

Given measurements of the probability density distri- 
bution, /, of a system, we solve for (or infer) the evolution 
equation of the system. Specifically, given the intensity 
profiles at zq, 2:2 and z&, we infer the complex Ginzburg- 
Landau equation 1(1)1. We start by considering the hy- 
drodynamic formulation of Eq. l|T)l. via the Madelung 




FIG. 1: The evolution of the complex field $ as a function of 
z. We use the intensity / at various slices Zo, Z2 and za, to 
infer the evolution equation of the system. 



transformation [lgj 



* = V7< 



;'<I> 



(4) 



where <3> = $(x,y) is the phase of the complex field VF 
Substituting Eq. J3J into Eq. (JTJ and, upon separating 
the real and imaginary components, we obtain two inde- 
pendent equations 



iv ± .(/v ± $) + ^|v ± $| 2 



G = 0, (5) 

- - 7 V X • UVx$ + - Vx® 2 + F = 0, 6 
az a 1 V / 2 



(2) where $ = 2$/a and 



1 az a a y/i 

2/(7) 2 1 _o /rp 

or a z y/I 



(7) 
(8) 



In this hydrodynamic form, Eq. © is a diffusion type 
equation with diffusion coefficient rj/a, whereas Eq. © 
is analogous to the Navier-Stokes equation for an com- 
pressible fluid [T^ . 

Inferring the TDCGL equation involves first solving 
Eq. 10 for the phase $ (see e.g., H3)- Because of the 
diffusive term, Eq. © becomes very complicated. Solv- 
ing for $ is thus nontrivial, even if 77, a and g(I) are 
known. 



III. PHASE RETRIEVAL IN THE PRESENCE 
OF DIFFUSION 

To illustrate the method in which 4> (or <E>) may be 
retrieved, we solve Eq. (J5J for $ assuming for the moment 
that a, 77 and g (I) are known. We write Eq. JSJ as 



Vl$ + V ± $ 



-G 



rja 



(9) 
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FIG. 2: (a) The exact phase at Z\ (A<s> = 0.1), (b) the phase 
retrieved from the simulated intensity data via the multigrid 
iterative phase retrieval scheme using the intensity at zo and 
Z2- White denotes a phase of 0.1 rad, whereas black denotes 
a phase of 0. This illustrates that it is possible to recover the 
phase of the system, given modulus information. 



The last term on the right hand side of Eq. © may 
be regarded as a perturbation term of second order in 
the phase gradient. This suggests a iterative numerical 
scheme to solve for the phase. At the kth iteration Eq. © 
is written as 



Vi$ fe + -V_l7- V ± $ fc 



-G 



V ± $ fc _i| 



(10) 



where is the phase at the previous, (k — l)th iter- 

ation. At the first iteration [k = 1) we set |Vj_$o| = 0. 
At successive iterations, is obtained using a multi- 
grid numerical scheme Il9], which is implemented by the 
Mudpack package p(Il2lll22l|. 

The iterative scheme assumes that the phase gradient 
is small. Therefore a successful and accurate recovery of 
the phase requires that |Vj_$| < 1. 

The convergence criterion of the iterative numerical 
scheme is determined by the change in the total phase 
gradient (the norm of the phase gradient) of the retrieved 
phase, ||Vj_$||, over two successive iterations. The 
scheme is said to converge when ||Vj_$||fc — ||Vj_$||fc-i 
is small, i.e., less than 10~ 6 . 

Figure|21shows a typical result; namely, the exact phase 
(the phase taken from the forward evolution) compared 
to the retrieved phase at zy. The exact phase at zy is 
obtained by averaging the exact phase at z and z 2 . At 
every iteration we monitor the Root Mean Square (RMS) 
error of the phase gradient, ct(|V^$|). This is defined 
according to 



<t(|Vx$| 



J n (\y±^\-\y^ r \fdxdy 

f n \V ± ^dxdy 



(11) 



where f2 € 3? 2 and $ r is the retrieved phase. Figure |3| 
(a) shows the RMS error for A<y, ~ 1, where the average 
phase gradient (|V X $|) = 0.982; Fig. |3| (b) shows the 
RMS error for A$ = 2, where the average phase gradi- 
ent (|V_l$|) = 1.97. The RMS error decreases rapidly 
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FIG. 3: The RMS error of the phase gradient, <r(|Vx$|), 
as a function of the iteration number k. (a) A$ — 1 (with 
(|V<E>|) = 0.982), showing that the error decreases quickly, 
and hence the numerical scheme converges, (b) A<s> = 2 (with 
(|V$|) = 1.97). Here the error oscillates illustrating that 
the numerical scheme does not converge. After about 390 
iterations the error in (b) diverges. 



in (a); however, the error in (b) oscillates and diverges. 
Numerical simulations show that the numerical scheme 
converges up to A<y> = 1.5 (i.e., (|V$|) ~ 1-47); how- 
ever, the scheme does not converge for A® > 1.6 (i.e., 
above (|V$|) ~ 1.57). This indicates that it is possible 
to retrieve the phase up to (|V<I>|) ~ 1.47. This is well 
above the constraint imposed by perturbative consider- 
ations, i.e., |V$| < 1. This suggests that the numerical 
iterative phase retrieval scheme is robust. The success- 
ful retrieval of the phase in the presence of diffusion is 
a significant step in our attempt to infer the complex 
Ginzburg-Landau equation. 



IV. PHASE SENSITIVITY WITH ERRORS IN 77 

The phase retrieval scheme is robust for a precisely 
known diffusion parameter rj. Here we investigate the 
effect of an error in rj on the retrieved phase. 

Figure 0] shows the comparison between the exact 
phase and the retrieved phase for various errors in rj. If 
we define the boundary of our system, to be a circle 
inscribed within the simulation domain [0,1] x [0,1], it 
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FIG. 4: Phase retrieval with various errors in the diffusion pa- 
rameter n. (a) Exact phase, (b) With no error in n, the phase 
can be retrieved accurately, (c) A negative error in rj gives 
the illusion of "particles" leaving the system, whereas (d) a 
positive error in n gives the illusion of "particles" entering the 
system. 



can be seen that the integral of the normal phase gradient 



TV 



Vi§ • ndL 



on 



(12) 



where dl is the path along the boundary, of the exact 
phase approximately vanishes on the boundary (see Fig. 0] 
(a)). The retrieved phase with an exact 77 as shown in 
Fig. 0] (b) is indistinguishable from the exact phase in 
(a). However, with negative error in 77, the quantity N is 
positive (see Fig. 01(c)). When the error in 77 is positive, 
N is negative (see Fig. 0] (d)). Since Vj_$ • n is the 
velocity of the "particles" entering or leaving the system, 
an N > in (c) gives the illusion of "particles" leaving 
the system; whereas • n < in (d) gives the illusion 
of "particles" entering the system. For the situation in 
which we want to infer the equation of evolution from 
modulus data, this suggests that accurate inference of 
the diffusion parameter 77 is required, if we are to infer 
an equation of evolution that accurately describes the 
system. 



V. INFERRING THE COMPLEX 
GINZBURG-LANDAU EQUATION 

In this section we describe how the TDCGL equa- 
tion Q may be solved for the unknown parameters a, 
77, and the unknown functions <E>(x,w), f(I) and g(I). In 
a manner similar to the way we infer the equation of evo- 
lution for the nonlinear dissipative Schrodinger equation 
[Io| . we divide the task into two consecutive parts. In 
the first part we solve for a, 77, $(x,y) and /(J) with 



known g(I). In the second part we discuss how g(I) may 
be determined. 



A. Inferring a, rj, <&(x,y) and /(/) 

The behaviour of N due to errors in 77 allows us to set 
up a "diffusion relaxation" iteration scheme to accurately 
infer a, 77, $(x,y) and /(/). We do this by assuming a 
priori knowledge of N; that is we know the current den- 
sity flowing through the boundary. This a priori knowl- 
edge of the boundary condition for the phase is not a se- 
rious limitation. For example, methods can be developed 
to measure the phase on the boundary. One such method 
was discussed in our previous paper (see PU). Here, for 
convenience, we assume N — on the boundary. This 
is satisfied, for example, if the system is trapped inside 
a potential well such as a Bose-Einstein condensate in a 
harmonic trap |23l |2 ll |25| . It is also satisfied if the sys- 
tem is confined to a container, such as for an uncharged 
superfluid [26|. In many other finite size systems, the in- 
tegrated current density normal to the boundary of the 
system is expected to be zero (if we choose a boundary 
that is sufficiently large). 

The diffusion relaxation scheme is as follows. In the 
first iteration we guess 77/a. This value is used to retrieve 
the phase $ and infer a. We then calculate N on the 
boundary. If N < 0, we know that the guessed 77/a 
is too large, otherwise it is too small. We then modify 
rj/a. The initial guess of 77/a can take any non-zero value; 
however, our initial guess is obtained by solving 77/a using 
Eq. (|SJ), with the assumption $>q(x,u) — 0. We solve for 
77/a (for the case &n(x, y) = 0) in a similar way to how 
we solved for a (see [l(|)- That is we substitute pairs 
of points with the same intensity, i.e., I\ — I2 (where 
Ii = I(xi,yi,z) and I2 = I(x2,y%, z)), into Eq. JSJ to 
obtain two independent equations - one equation for each 
of the points {x\,yi) and (£2, "2)- Subtract these two 
equations from one another to eliminate g(I) and obtain 



V 
a 



dh 



dl 2 

dz 



(13) 



By using many constant-intensity surfaces a histogram 
of 77/a can be constructed, and 77/a is obtained as the 
peak of the histogram (in the specific case when g (J) is 
known, Eq. JSJ may be solved directly for 77/a). 

To guide us regarding how to modify 77/a at successive 
iterations, we can guess 77/a again at the second iteration. 
For example, at the second iteration we increase 77/a by 
1% of the initial guess. This results in N in the second 
iteration differing slightly from that at the first iteration 
by a fractional amount X = (772 — f}i)/f}x (where the 
subscripts denote the iteration number). The magnitude 
and direction in which N varies allows us to systemat- 
ically modify 77/a at the next iteration. This diffusion 
relaxation algorithm is summarised in Appendix IA II 

The diffusion relaxation method fits well with the iter- 
ative phase retrieval method that is used to retrieve the 
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phase of the complex system. For example, at the first 
iteration a is unknown (even though 77/a was approx- 
imated or guessed initially). However, it is immaterial 
since initially we set <&o(x,y) = and consequently the 
term ^-|V^<I| 2 in Eq. JSJ vanishes. At successive iter- 
ations 2 ^|Vj_ < I > | 2 is approximated by using $ from the 
previous iteration. 

In the diffusion relaxation scheme, N at two succes- 
sive iterations is used to modify 77/a at the next iter- 
ation. The algorithm to find the fractional increase, 
X = (r?fc+i - Vk)/Vk, in y/a is given in Appendix IA 21 
The behaviour of the algorithm is as follows. If ./V changes 
from negative to positive (or vice versa) at two succes- 
sive iterations, X changes sign. In this situation, the 
magnitude of X is smaller than its previous value. This 
guarantees that the numerical scheme is stable and that 
N evolves towards zero. If N at two iterations has the 
same sign, X negative (positive) (i.e., 77/a will decrease 
(increase)) depending on whether N deviates away from 
zero (or approaches zero). The algorithm for obtaining 
N is analogous to the Newton-Raphson method for find- 
ing the root of a nonlinear function |27j . Consequently, 
the convergence of this algorithm is robust. At the end 
of each iteration, we use X to update rj/a. 

For each iteration \X\ is monitored to test for conver- 
gence. If \X\ < e, where e is a small tolerance, we con- 
sider the numerical scheme to have converged. Since the 
retrieved phase is sensitive to errors in 77 in the order of 
~ 0.01, the tolerance e is set such that rjk+i ~ Vk <C 0.01, 
i.e., e = 10~ 7 . During the iteration scheme 77/a evolves 
towards an asymptotic value, (77/0)00. 

To ensure that the numerical scheme converges from 
both side of (77/01)00, we reflect the initial guess of rj/a 
through (77/0)00 and start the iteration scheme again 
(with $(a;,y) reset to zero). The final diffusion param- 
eter, 77, is taken as the average of that obtained from 
the two iteration schemes. This average value is used to 
retrieve <E>, compute o and to retrieve the function /(/) 
using Eq. ©. 

Figure 13(a) shows the evolution of 77 as a function of 
the iteration number k. For an arbitrary guess of 77/a, the 
numerical iteration scheme evolves 77 towards the asymp- 
totic value 77 — 2.00044. Furthermore, when the initial 
guess is reflected through this asymptotic value and the 
iteration scheme re-run, 77 tends again towards the same 
asymptotic value. This illustrates that the numerical 
scheme is robust and 77 converges quickly. The actual 
value of 77 is 2. The error in 77 may be due to the integral 
of the normal phase gradient not precisely vanishing on 
the boundary and from the error in the retrieved phase. 

Figure (b) shows that o also evolves towards an 
asymptotic value, o = 1.007. When 77 is overestimated, 
the parameter o evolves towards the asymptotic value 
from below; however, o is less sensitive when 77 is under- 
estimated. The actual value of o is unity. The error in 
o may also be due to the error in 77 and the error in the 
retrieved phase. 

Figure (c) shows the behaviour of N for N > and 
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FIG. 5: The evolution of 77, a, N and X in the diffusion 
relaxation numerical iteration scheme for a typical simula- 
tion (d'l'j} — 0.49). (a) 77 converges to a precise value in the 
long term evolution. From above 77+ — * 2.000443, from below 
?7_ — > 2.000439. The exact value is 77 = 2. (b) When 77 is over- 
estimated, a was initially underestimated (its exact value is 
unity); however it subsequently evolves towards q_ — > 1.007. 
When 77 is underestimated, a is relatively constant illustrating 
that a is less sensitive to negative error in 77, when compared 
to a positive error in this quantity, (c) shows that the normal 
phase gradient along the boundary tends towards zero from 
both sides, whereas (d) illustrates the fractional variation of 
77 as it converges towards a solution. 



for N < 0. The magnitude of N is large at the start 
of each iteration. Its values tend towards zero quickly 
from both sides as a function of iteration number, illus- 
trating the convergence of the numerical scheme. The 
corresponding fractional change in the diffusion parame- 
ter, X , plotted in Fig. 0(d) indicates how 77 varies during 
the iterations. The parameter X is varied such that N 
evolves towards zero, and 77 and o evolve towards their 
asymptotic values. 

A typical frequency histogram of o at the end of the 
iteration is shown in Fig.0((|<I>|) = 0.49). The value of o 
taken from the peak of this histogram is o = 1.01 ± 0.04, 
with the error taken as the Full-Width-at-Half-Maximum 
(FWHM) of the frequency histogram. The plot of the ex- 
act phase and the retrieved phase at the end of the sim- 
ulation are shown in Figs. [7] (a) and (b). No discernible 
difference is found between Figs. [7] (a) and (b). The frac- 
tional RMS errors in the phase and in the phase gradient, 
cr(<i>) and er(|V^<i>|), are shown in Figs. 0(c) and (d). In 
both cases the fractional RMS errors are less than 2%. 
Note that even when we use the exact 77 and o, the frac- 
tional RMS errors in the phase and in the phase gradient 
are also of a similar value. 

Once 77, o and $ are obtained, we can calculate the 
nonlinear term /(/). This is a straightforward process, 
which involves the direct application of Eq. JSJ . A typical 
result is shown in Fig. Besides a constant energy shift 
of —131.14 (see 01), ^he error in the inferred nonlinear 
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FIG. 6: A typical frequency histogram of a at the end of the 
diffusion relaxation scheme, for which a = 1.01 ± 0.04. The 
error is taken as the FWHM of the histogram. The histogram 
is constructed from 100 equally spaced iso-intensity surfaces 
using Eq. ©. The vertical axis is the number of occurrences 
of a within the interval Aa = 10 -3 , whereas the horizontal 
axis shows the various retrieved values of a. 



FIG. 8: Comparison between the exact function /(/) = 
100sin(7rJ) (upper curve) and the inferred function (lower 
curve), which can be fitted with —131.14 + 100.59 sin(7r/). 
This illustrates that the nonlinear function /(/) can be re- 
covered accurately. The constant off-set is indicative of the 
phase, which can only be retrieved up to an arbitrary con- 
stant. 
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FIG. 7: Comparison between the exact phase (a) and the 
retrieved phase (b). (a) and (b) show no discernible difference 
between the exact phase and the retrieved phase at the end 
of the diffusion relaxation scheme, (c) and (d) show that the 
RMS error in the phase and the phase gradient quickly decays 
with increasing iteration number k. 



term /(/) is small and can be accounted for from the error 
in inferring r\ and a, and from the error in the retrieved 
phase <5. 

So far we have neglected the dissipation term g(I). The 
diffusion relaxation scheme assumes that the functional 
form of this dissipation is known. In Sec. IV Bl we discuss 
the effect of dissipation on the phase of the complex field 
'J, together with ways in which the nonlinear dissipation 
function g(I) can be inferred. 



B. Inferring g(I) 

Our diffusion relaxation scheme assumes that the dis- 
sipation function g(I) is known. In this section we show 
that, in general, it is hard to solve for the dissipation 
in the same numerical scheme that is used to infer the 
parameters rj and a, and the functions /(/) and $ (ex- 
cept in the case where &(x,y, z) — 0). This is because 
dissipation has an adverse effect on the evolution of the 
phase, in the sense that will now be described. 

Consider the case when | V$| is small, for which we can 
write Eq. JSJ as 



dz 



V 



_L • 



(/Vil) 



2Ig(I) 2 V 



_^V7V 2 V7«0. (14) 



If we define 



V • (IV H) EE -g(I), 
a 



(15) 



where H = H(x, y, z) is a real-valued function of posi- 
tion, Eq. I|14l) may be written as 



— + V 

dz 



IV j_ $ + H 



-—VlV 2 Vl~0. (16) 
a 



Equation (| 1 fc>|) is a continuity equation in the presence of 
a new phase distribution $ + H. In this equation, the 
dissipation term completely disappears. Since Eq. l(TB|) 
is identical to Eq. (|14|) . we see that dissipation modifies 
the "phase" of the wavefunction from $ to $ + H, or 
equivalently changes the "flow" velocity by V±H. The 
new current density is /Vj_($ + H). This implies that 
dissipation is closely related to the phase of the complex 
field and one cannot distinguish between the two entities, 
using the scheme previously outlined. 
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The phase <E> can only be determined if H is known or 
vice versa. Therefore to measure all terms in the TDCGL 
equation, the dissipation of the system either has to be 
measured separately (i.e., separately prepare the system 
so that the dissipation can be measured without a priori 
knowledge of other parameters and functions), or it has 
to be measured by other means. In \ld^ we discussed 
two ways in which the dissipation function g(I) can be 
measured. Here we generalize to include diffusion. 

In many systems, such as monoenergetic electron 
beams or electromagnetic waves, it may be relatively easy 
to prepare a plane wave state with a constant transverse 
intensity profile in which the phase <E> is independent of 
the x and y positions, over the plane at constant z. Sup- 
pose we can prepare the system in such a plane wave 
state. In this situation the divergence of the current den- 
sity Vj_ ■ (7V_l$), the second order velocity field |Vj_$| 2 
and the diffraction term I~ 2 V 2 -\/7 vanish. Using Eq. JSJ, 
the dissipation is then given by 



9(1) 



a 



1 dl 
"2Idz' 



(17) 



Equation (|17|) is independent of the diffusion parameter 
7]. This illustrates that for a uniform plane wave, the 
system does not diffuse. To obtain the dissipation at 
different values of the intensity, we need to make repeated 
measurements. We can separately prepare the system 
for each of the intensity values, or since the intensity 
always decreases due to dissipation, we can measure the 
dissipation over long time evolution. For the latter case, 
we note that the intensity evolves as 



I(z) = I e-i $ 9{z)dz , 



(18) 



where Jo is the intensity at z — 0, indicating that dissipa- 
tion leads to decay in the intensity profile of the system. 

In some systems, such as an uncharged superfluid or a 
Bose-Einstein condensate, it may be difficult to construct 
a plane wave state. In such a case we examine an alter- 
native approach to measuring dissipation. As outlined 
in [ljj for the non-diffusive case, an alternative approach 
is by averaging the measured dissipation over sufficiently 
large number of measurements, M. For the diffusive case 
discussed here, we restrict ourselves to systems with low 
fluctuations so that |Vj_$| is small. In such cases, the 
second order term in the velocity field, |Vj^$| 2 , is neg- 
ligible. The continuity equation is then reduced to that 
given by Eq. i|14|) . The method of averaging discussed 
in ^(j is applicable and the dissipation can be approxi- 
mated as 



a 



M 

-Y 



k=l 



a 1 vV7 _i^ 

a yfl 21 dz 



(19) 



where k denotes the kth measurement. In our discus- 
sion here, r\ja is not known. So the dissipation func- 
tion should be measured in conjunction with measuring 
all other parameters and functions in the TDCGL equa- 
tion (Q. That is, the dissipation should be measured in 



the same iteration scheme that is used to measure a, 77, 
<& and /(/). The dissipation function may be calculated 
immediately using Eq. 1)19(1 after rj/a has been found or 
approximated at each iteration. 



VI. GENERALIZATION TO 
MULTI-COMPONENT COMPLEX FIELDS 

The diffusion relaxation algorithm for inferring the 
TDCGL equation may be generalized to the case of multi- 
component (2+ l)-dimcnsional complex fields, denoted by 
{& n (x,y, z)}, which comprise a set of N complex fields 
ee V n (x,y,z),n = 1,---,N. The TDCGL equation 
governing the evolution of this multi-component complex 
field may be written as 



Q 

ia nTT + (1 - *?/n)V 2 L + /„ 
OZ 



= 0, (20) 



where a n and rj n are real numbers, while f n (I\,' ■ ■ , In) 
and g n (h, ■ • ■ > In) ar e real- valued functionals dependent 
on the intensity /„ = l^nl 2 - The "hydrodynamic" for- 
mulation of Eq. (|2*U|l is: 



In 



dz 



|V ± $ n | 2 + G r , 



0, (21) 



— ^Vj.- (l„V x $„ 

+ ^|V ± $„| 2 +K = 0, (22) 



where $ n = 
G„ 

Fn 



2arg(* n )/a n and 

1 dl n 2gn _ 

I n dz a„ 



Vn 

2 



a 



fin 



2r? n 1 

' fT n . 



(23) 
(24) 



Given modulus information on three closely spaced 
planes, Eq. (|21|l may be used to solve for r]n/a n and 
<!>„, using the diffusion relaxation iteration scheme, sepa- 
rately for each component. However, to solve Eq. (|22ll for 
a n , our technique is to find pairs of points with the same 
f n . For a single component complex field, such pairs of 
points can be found on iso-intensity surfaces; however, 
this is not always the case for multi-component complex 
fields. For arbitrary multi-component complex fields it 
is not known how such pairs of points can be found (al- 
though such pairs of points exist since /„ vanishes on the 
boundar y an d is non- vanishing in the interior). We out- 
lined in |lOj that, in principle, for two-component com- 
plex fields (N = 2) in two spatial dimensions, we can 
always find such pairs of points. Our finding in [ljj is 
also applicable to the TDCGL equation considered in 
this paper. Here we generalize the arguments for find- 
ing pairs of points with the same f n for two-component 
complex fields (N = 2) in two spatial dimensions. We 
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y 




x Ti 
(a) (b) 

FIG. 9: Schematic illustration of the trajectory Ti and Tb 
corresponding to a path in Ji and /2-space. (a) Ti traverses a 
path in physical space corresponding to an iso-intensity sur- 
face in Ji. (b) Along the same trajectory in physical space, 
the intensity I2 varies in T2, with the end points identified. 
In general, the trajectory T2 is not an iso-intensity surface; 
however, there are pairs of points with the same intensity. 



cal model is required to explain the physical phenomenon 
and make sense of the experimental data. 

Knowledge of the evolution equation of a system not 
only allows us to obtain quantitative understanding and 
physical insight into the system, but also to study the 
future evolution of systems where long term observations 
are not possible. Future investigations will be directed 
towards extending the TDCGL model to include gauge 
fields, which will allow us to describe larger classes of 
systems including those which exhibit superconductivity. 
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construct a closed trajectory, T n (n = 1 say), in which 
I\ is an iso-intensity surface (see Fig. 0(a)). As we tra- 
verse a path in T\, 1% traverses the corresponding path 
in T 2 as shown in Fig. [3] (b), where not every point on 
T 2 necessarily has the same intensity. However, since T\ 
is a closed trajectory, T 2 is necessarily a closed trajec- 
tory, i.e., the two end points of trajectory T 2 in Fig. |5| 
(b) are identified. Consequently there are at least two 
points in T 2 with the same intensity, and it is easy to 
follow the methods developed in this paper to infer the 
TDCGL equation of the two-component field. This ar- 
gument can be generalized to infer the TDCGL equation 
for a three-component complex field in three spatial di- 
mensions. In three spatial dimensions, it is possible to 
construct a closed iso-intensity surface T\ for 1\, For 
some paths in Ti, it is possible to trace out a closed iso- 
intensity loop T 2 for I2, and for some points in T 2 , it is 
possible to find pairs of points with the same intensity 
for I3. Therefore in three spatial dimensions, it is possi- 
ble to infer the TDCGL equation for a three-component 
complex field. 



VII. CONCLUSION 

We have presented a novel method for 'measuring' the 
evolution equation of a two-dimensional complex wave 
field obeying the TDCGL equation, given only the mod- 
ulus information of the wave field. This is done via a 
diffusion relaxation and a Newton-Raphson type itera- 
tive scheme. The numerical scheme is robust; however, 
it is only applicable for high signal to noise ratio in the 
modulus information, i.e., higher than 10000:1 for our nu- 
merical resolution. This presently restricts the applica- 
bility of our methodology to systems where it is possible 
to make very precise measurements on the system. 

Notwithstanding the noise problem, this work may be 
of general significance in a variety of fields, where a physi- 



APPENDIX A: NUMERICAL METHODS 

Here, we give details of the numerical scheme used to 
infer the TDCGL equation based on the diffusion 
relaxation scheme discussed in Sec. IV Al The diffusion 
relaxation scheme is given in IA II Within the diffusion 
relaxation scheme a Newton-Raphson type convergence 
algorithm, given in I A 21 has been implemented to modify 
the diffusion parameter at each iteration. See main text 
for the definition of symbols used. 

1. Diffusion relaxation scheme 

•Approximate (or guess) r\ja 

•Set iteration number k = 1 

•Retrieve $ at Z\ and z 3 using Eq. (J5J 

•Solve for a on iso-intensity surfaces using Eq. (j^J) 

•Calculate ATj. (N at the kth iteration) 

•Increase rj/a by, say 1% 

1. set k = k + 1 

2. retrieve <E> at z\ and 23 using Eq. J5| 

3. solve for a on iso-intensity surfaces using Eq. JSJ) 

4. calculate N k 

5. use (i]/a)k, Nk and Nk-i to find a new r]/a 

6. if \Ar)\/r] > e (a small tolerance), repeat step 1. 

2. Diffusion modification algorithm 

IF N k+1 N k < THEN 

X ~~N k + 1 -N k X 

ELSE 

IF |V fe+1 | > |JV fc | THEN 



N k+1 
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END IF 
END IF 

7] k+1 = (1 + X) Vk 
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